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Abstract. Nowadays sp carbon chains terminated by graphene or graphitic-like 
carbon are synthesized routinely in several nanotech labs. We propose an ab-initio 
study of such carbon-only materials, by computing their structure and stability, as well 
as their electronic, vibrational and magnetic properties. We adopt a fair compromise 
of microscopic realism with a certain level of idealization in the model configurations, 
and predict a number of properties susceptible to comparison with experiment. 

PACS numbers: 81.07.-b, 31.15.A-, 75.70.Ak, 65.80. Ck 
1. Introduction 

The carbon atom, with its three possible hybridization states, originates in nature very 
different elemental materials. The three possibilities (sp 3 , sp 2 , sp) correspond to three 
different prototypical structures: respectively diamond, graphitic-like structures (such as 
graphite, graphene, carbon nanotubes, and fullerenes), and linear carbon chains (known 
in the literature as polyynes [U E] or sp carbon chains, spCCs in short). 

spCCs were discovered in nature around 1968 [3J, but their role in the arena 
of carbon-based nanostructures has been quite marginal: till very recently, spCCs 
have been considered as exotic allotropic forms, mainly present in extraterrestrial 
environments. Indeed, in the interstellar clouds formed in the explosions of carbon stars, 
novae, and supernovae, spCCs have been detected aside with fullerenes (HE], amorphous 
carbon dust, cyanopolyynes, and oligopolyynes. Indeed, in the phase diagram of carbon, 
the field of existence of spCCs coincides with that of fullerenes [6] . 

The high reactivity of spCCs [7], and their tendency to undergo cross-linking to 
form sp 2 structures, directed the experimental efforts towards complicate strategies 
for the stabilization of the spCCs with molecular end-groups or their isolation in 
inert matrices [D El E|. Thanks to novel synthetic routes and strategies [TU], 
polyynes of increasing length and different type of termination have been successfully 
synthesized and characterized [TH [121 EES EH ES]. End-capped spCCs (often, but not 
exclusively, hydrogen-terminated) are being synthesized by chemical / electrochemical / 
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photochemical methods [U HE! [TTJ [TBI HE] . spCCs have also been produced from carbon 
by dynamic pressure (TTl [20] . spCCs or carbynoid material containing up to 300 carbon 
atoms were synthesized, which opens a promising route toward molecular engineering 
of sp-carbon structures [HI [TTJ [211 122] . 

Samples of pure carbon films grown by supersonic cluster beam deposition at room 
temperature have been characterized, and proven to contain a sizeable sp component 
[23] 121] . These spCCs showed weaker stability relative to the sp 2 component upon 
exposition to low pressure gases at room temperature [25 J . spCCs have also been 
produced from stretched nanotubes [26j. Recently, Jin et al. [21] have realized spCCs 
by stretching and thinning a graphene nanoribbon from its two free edges, by removing 
carbon rows until the number of rows becomes one or two. These spCCs show a good 
stability under the beam of a transmission electron microscope (TEM) for lengths up 
to a few nanometers. These experimental observations of spCCs formation during 
the controlled electron irradiation of graphene planes resulted in an rapidly increasing 
interest in this field [2SII2SIEQIEIII321E31E11I351I35]. 

Meanwhile, the existence of intrinsic magnetism in pure carbon has been a matter 
of debate for quite some time [37]. Possible effects due to magnetic contaminants on 
the experimental results have been discussed [38]. Subsequently, it has been shown 
that possible contamination effects are unable to explain quantitatively the measured 
ferromagnetism, supporting the idea that carbon magnetism has an intrinsic origin [39J . 
The existence of 7r-electron magnetism in pure carbon has now been widely accepted 
(see e.g. [101 Hi]). On the theoretical side, it is well known that magnetic instabilities 
exist at specific graphene edges [121 EH EH H5] , in defective graphene [46J and nanotubes 

The discovery of pure-carbon 7r-electron magnetism has also lead to speculations 
about possible applications of carbon-based magnetic materials in molecular electronic 
devices: for example, spintronic devices built around the phenomenon of spin- 
polarization localized at the 1-D zig-zag edges of graphene have been proposed [48J. 
Recently, the accent has been also put on the possibility to modify the magnetization 
optically [19]. Indeed, linear carbon chains are nowadays considered promising 
structures for nano-electronic applications [SOI EI]- F° r example, they can be used 
as molecular bridges across graphene nanogap devices. Potential applications for the 
realization of non-volatile memories and two-terminal atomic-scale switches have been 
demonstrated [5l] . 

The remarkable robustness of spCCs terminated on pure carbon, graphene-like 
fragments, combined with the fast progresses in the synthesis of graphene and graphene 
derivatives, could open the way towards the realization of actual nanodevices based on 
sp + sp 2 carbon nanostructures. The interest of such systems stems from the possibility 
to exploit their peculiar semiconducting-magnetic behavior, in order to achieve novel 
characteristics and functions for the target devices. 

Clearly, the possibility of designing graphene-based magnetic nanostructures is 
particularly intriguing. The capability of arranging the spins inside a carbon structure 
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in a variety of ways, could open the way for the construction of completely novel devices 
[52J. Possible future applications for spCCs in interaction with the graphene-type 
system could be the construction of microchips with ferromagnetic or antiferromagnetic 
character that can be controlled by nanomanipulation and read out by nanocurrents. 

In Sect. [2]we introduce the theoretical model, whose electronic structure we address 
by a standard ab-initio method based on the Density- Functional Theory (DFT). Section 
[3] presents the investigation of the structural and binding properties of a spCC inside 
a nanohole (nh) in a graphene sheet representing the sp 2 component in a carbon- 
only sample. Section [| studies the magnetic properties of the nh edges and of the 
inserted spCCs. In Sections [5] and [6] we cover the band structures and selected 
vibrational properties of the studied nanostructures. In Sect. [7] we investigate the 
dynamical stability of the metastable spCC-nh structures, by means of a tight-binding 
molecular dynamics model. Section [8] discusses the results of simulations in the light of 
experimental data. 

2. The model 

In the present paper, we focus on spCCs bound to graphitic structures, represented by 
a hole in a infinite graphene sheet. This system is representative of a class of sp + sp 2 
systems, which are at the core of an intense experimental work (231 I2H [221 123 12E1 EIH EH 
155] . To address the structural, vibrational, and electronic properties of spCCs inserted 
into a nanometer-sized hole defect of a graphene layer, we resort to the DFT. The plane- 
wave pseudopotential method and the Local Spin-Density Approximation (LSDA) to 
DFT have provided a simple framework whose accuracy have been demonstrated in a 
variety of systems |56j. The time-honored LDA is one in many functionals being used 
for current DFT studies of molecular and solid-state systems: other functionals often 
improve one or another of the systematic defects of LDA (underestimation of the energy 
gap, small overbinding and tiny overestimation of the vibrational frequencies), but to 
date no functional is universally accepted to provide systematically better accuracy than 
LDA for all properties of arbitrary systems. For a covalent system of s and p electrons as 
the one studied here, LDA is appropriate, and we expect our results to change by a few 
percent at most if the calculations were repeated using some other popular functional 

[HUES EH]. 

We compute the total adiabatic energies by means of the code Quantum Espresso 
[60J, which computes forces by standard Hellmann-Feynman method. Each self- 
consistent electronic-structure calculation stops when the total energy changes by less 
than 10~ 8 Ry. We use ultrasoft pseudopotentials [M1E2], for which a moderate cutoff 
for the wave function/charge density of 30/240 Ry is sufficient. We terminate atomic 
relaxation when all residual force components are smaller than 10 -4 Ry/a ~ 0.04 /iN. 

Plane waves require periodic boundary conditions in all space directions. In our 
model, we represent a graphene plane in a (x — ^-periodically repeated 1714 x 1690 pm 2 
supercell consisting of 7 x 4 rectangular conventional unit cells containing four carbon 




Figure 1. (Color online) The adopted slab model supercell for graphene, composed 
by 7 x 4 conventional rectangular cells containing 4-atoms, each equivalent to two 
primitive unit cells of graphene (in gray). Solid: the unit vectors of the conventional 
cell, of length a — 244 pm and b = ay/3 ~ 422 pm. The size of the full supercell 
(delimited by the solid lines) in the x — y plane is 1714 pm x 1690 pm. 

atoms each, see Fig. [Tj Within this cell, a complete graphene layer is then represented 
by 112 carbon atoms. We represent the graphite surface in the slab approximation, as 
one or a few stacked graphene layers: to ensure that the interaction between periodic 
images of the graphene sheet is negligible we interpose at least 1 nm of vacuum along 
the z direction. We have selected this cell size as it allows us to create reasonably large 
nanoholes in the graphene layers, with fairly small interaction between the supercell 
repeated copies of the nanohole itself, at the price of a manageable computer time. We 
cover electronic band dispersion the electron bands in the horizontal plane by means of 
a 5 x 5 x 1 T-centered k-point mesh. BZ integration of the metallic band energies are 
performed using a 2 x 10 -4 Ry-wide Gaussian smearing of the fermionic occupations. 

2.1. The Nanohole 

Starting from the perfect graphene foil of Fig. [TJ we remove selected atoms in order to 
form a nh. The size of the nh should be such that inserted spCCs fit and only bind at 
their ends. If the nh is too small, then spCC atoms would tend to reconstruct the edges 
of the hole, as illustrated in the examples of Figs. [2] and |3j There, the small spCC-edge 
distance leads to spontaneous (barrier-free) edge reconstruction with the formation of 
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(a) Initial (b) Relaxed 



Figure 2. (Color online) Initial and relaxed positions of a small nh with a G@ polyyne: 
the lateral nh size is so small that no barrier prevents two polyyne atoms (brown/clear) 
to bind to the armchair edges. 




(a) Initial (b) Relaxed 



Figure 3. (Color online) Initial and relaxed positions of a small nh with a Cq polyyne: 
the lateral nh size is so small that no barrier prevents two polyyne atoms (brown/clear) 
to bind to the zig-zag edges and form pentagonal rings. 
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Figure 4. (Color online) A nh of sufficient size to lodge spCCs is obtained by 
removing 28 atoms from the perfect graphene of Fig. [T] The rectangular hole has size 
975 pm x 985 pm, and is delimited by both zig-zag and armchair edges 



additional squares or pentagons. 

We must therefore construct a large enough nh, in particular allowing for a distance 
of at least 300 pm between the spCC and the nearest nh edge in order to prevent 
recombination reactions. Figure [4] shows the minimal nh with such property. Starting 
from the perfect graphene of Fig. [TJ the nh is obtained by removing 28 atoms forming 
a rectangular hole of size 975 pm x 985 pm with both zig-zag armchair edges. The size 
of the nh permits us to investigate the insertion of several spCCs, from C5 to C§ atoms 
long, with various positions relative to the nh edges. 



3. spCCs Binding to a Nanohole 

In an experiment [23J, |2U I2S1 E2] , spCCs are likely to be bound to extended sp 2 structures, 
more reminiscent of graphite than of graphene. Such a configuration could be described 
for example with one layer of perfect graphene plus a layer with a nh, such as that 
described above Sect. |2.1| Figure [5] shows this configuration with the insertion of a Cq 



polyyne. The upper graphene layer and the polyyne are fully relaxed, while the lower 
layer is kept frozen in ideal graphitic positions. The equilibrium distance between the 
two layers, 332 pm, is very close to the observed interlayer distance of graphite, 335 pm. 

The resulting system composed of over 200 atoms is computationally quite 
expensive: indeed, with a modern parallel computer using 32 Xeon-class processors, 
it took more than two weeks to obtain a fully relaxed configuration. On the other 
hand, we verified that the properties of the upper layer do not change significantly 
without the lower (perfect) sheet, whose effective corrugation is quite small. Indeed, 
the forces between the two layers are weak long-range forces whose action is mild and 
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Figure 5. (Color online) A nano-indented graphite surface represented by a perfect 
graphene foil (pale gray) plus a superposed sheet with a nh (dark blue) as in Fig. [4] 
stacked in the standard AB arrangement of graphite. A Cg chain (also dark blue) is 
connected to the zig-zag edge of the nh. The layer-layer distance equals 332 pm. 

almost translationally invariant, compared to the intra-layer forces acting in the nh 
layer. The DFT-LSDA evaluation of such weak dispersive forces is unreliable anyway. 
These observations suggest that it makes sense to consider the single layer containing 
the nh and the C n spCC inserted into it, and leave the fixed substrate layer out. The 
relaxation of the positions, performed in the same conditions as above, took about one 
week only in this 1-layer configuration involving 90 C atoms rather than 200 ones. 

In our calculations, we consider several spCCs, from C 5 to C 8 , placed in different 
positions inside the nh. We identify such compounds as nh-C n , with further specification 
when different local minima are considered. Their relaxed configurations are depicted 
in Fig. [6} Table [T] summarizes the structural properties of the configurations considered, 
comparing in particular their bond-length alternation (BLA) [63] and the spCC- 
graphene bonding energy, defined as the total energy of the empty nh plus that of 
the isolated spCC minus the total energy of the bonded spCC-nh configuration under 
consideration. 

For the selected nh size, C n chains of different length can fit more or less easily 
inside the nh. Short chains such as C5 or may fit at the price of a tensile stress, 
while longer spCCs can be forced inside the nh with a compressive strain, which could 
be eased by buckling [3U EU [65] . 

Strain influences directly the spCC BLA: a tensile strain leads to stretching more 
the weaker bonds, thus producing an enhanced BLA, typical of polyynic spCC (in nh-C6 
arm, whose BLA reaches 12 pm). Likewise, a BLA~ 10 pm is obtained for nh-Cg zig 
due to the nh size being approximately 5% larger than the equilibrium length of the 
spCC. In contrast, a compressive strain leads to a more cumulenic-type structure, e.g. 
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(g) nh-C7 curved, side view (h) nh-C7 s-curved, side view 




(i) nh-C7 straight, side view (j) nh-Cs side view 




0.0 0.5 1.0 1.5 nm 

Figure 6. (Color online) The relaxed configurations of the spCC-nanohole structures 
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[Fig-] 



Description 



BLA 
[pm] 



Ehond 

[eV] 



<7 r , 



Oh 



nh-C F 



[Fig. 6(a) 



nh-C fi lb 



[Fig. 6(b) 



A C5 chain stretched across the nh bonding 
weakly to both zig-zag edges (bond lengths: 
170 pm). 

The C 5 chain bonded to one zig-zag edge of 
the nh; all bond lengths of the spCC are 
similar to the lengths of cumulenic double 
bonds (~ 127 pm). 
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4.2 



6.2 



nh-C6 zig 



[Fig. 6(c 



A weakly stretched Cq chain bonded to the 
zig-zag nh edge. 
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12.9 



nh-Cfi arm 



[Fig. 6(d) 



The Cg chain joining opposite armchair nh 
edges. 



12 



8.2 



nh-C7 curved 



[Fig. 6(g 



nh-C? s-curved 



[Fig. 6(h) 



nh-C7 straight 



[Figs. 6(e) , 6(i 



A C7 connected to the zig-zag edges, and 
buckling out of the graphene plane (maximum 
spCC height ~ 120 pm). 
A C7 connected to the zig-zag edges, buckling 
in a s shape, with the central atom in the same 
plane as graphene. 

A compressed straight C7 joining the zig-zag 
edges. 



12.0 



11.9 



11.8 



nh-C, 



[Figs. 6(f). 6Q) 



A Cs curved chain. The maximum height of 
the spCC equals 297 pm. 



12.5 



wnh-2Cg 
[Fig. [7] 



Two Cfi spCCs inserted in a wider nh joining 
the zig-zag edges. The lateral distance 
between the spCCs equals 491 pm. 
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n.c. 



Table 1. (Color online) Summary of the individual configurations considered for the 
spCCs bound to the nh. For the relaxed configurations we report the resulting BLA 
and total bonding energy Et, on< i, corresponding to the formation of the (usually two) 
bonds between the spCC and the nh. The relevant symmetry planes (<r p is the reflection 
across the graphene layer plane, a v is the vertical plane through the spCC, cr^ is the 
horizontal plane through the middle of the spCC) are marked for the configurations 
for which they apply. 

*Individual a p and ah are not symmetries for nh-C7 s-curved, but their product a p <Th 
is. 
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Figure 7. (Color online) Relaxed position of two Cq spCCs in the wider nh: the hole 
width is 1504 pm, large enough to prevent the recombination of the two spCCs. This 
configuration displays no significant novelty relative to the single polyyne nh-Cg zig 
configuration. 

1111-C7 straight has BLA ~ 2 pm. As was observed in a slightly different context |34j, 
very small BLA variations are induced by lateral atomic displacements. 

The linear size of the hole is approximately 95% of the equilibrium length of 
the C7 chain: different stable shapes of the spCC in nh-C7 can be stabilized by the 
compressive strain [64]. We study three equilibrium geometries of the spCC: straight, 
single-curvature buckling, and s-curved buckling. The compressive strain depresses the 
BLA, so that for the three of them the BLA ranges from 2 to 3 pm. 

The nh-Cs is so much stretched that if kept in a central symmetric configuration 
bonding between the spCC and the two edges of the nh is weak, each highly stretched 
terminal bond contributing only about ~ —2 eV to lowering the total energy. In such a 
condition we observe an intermediate BLA ~ 7 pm. This configuration is locally stable, 
but if we displace the spCC significantly (~ 50 pm) closer to one nh edge than to the 
other, and then let it relax, we retrieve an energetically favored configuration (nh-Cs 
16) with essentially a single strong bond (total energy lowering: ~ —6 eV) between the 
spCC and the nh. Here the spCC internal bond lengths are practically equal to those 
of isolated C5. 

Due to the small size of the nh, the C§ chains can only fit in a curved geometry: the 
maximum out-of-plane elevation of the spCC equals 297 pm. The resulting BLA ~ 6 pm 
is intermediate between cumulenic and polyynic. 

We also consider a wider nh in which one can insert more than one spCC: in wnh- 
2C6, the nanohole contains two C6 spCCs at a distance large enough to keep them 
separated, see Fig. [7j All properties are essentially equivalent to those of the nh-C 6 zig, 
therefore we will not further investigate this configuration. 

We evaluate the bonding energy of the configurations described here. Due to 
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Figure 8. (Color online) Scheme of the geometric relation between adjacent graphene 
edges. The angle between edges is defined as the angle between the vectors normal to 
the edge. The edge atoms belong to the same sublattice (either dark/red or clear/gray) 
when the zig-zag edges are at a relative angle of 0° or 120°; they instead belong to 
different sublattice when the relative angle is 60° or 180°. (Adapted from Ref. [57].) 



its stretching, the nh-Cs has a little value of Ebond ~ 4 eV, while for the 11I1-C5 lb 
Ebond = 6.2 eV which can be considered a fair estimate of the spCC-graphene edge 
binding energy according to DFT-LSDA, and matches previous evaluations [25]. For 
all other configurations Eb on d — 12 eV indicative of the formation of two bonds, at the 
expense of approximately 1 eV which accounts for the elastic deformation energy of the 
spCC and the connected graphene. 

Until now the spCC was always connected to zig-zag edges. When a C 6 chain binds 
to the armchair edges (nh-C6 arm), the bonding energy is smaller (Ebond — 9 eV), due 
to the lower reactivity of the armchair edge relative to the zig-zag one [25| 166] . The 
BLA assumes a highly dimerized value 12 pm associated to a tensile strain, like for the 
nh-Cg zig isomer. 

In the following we shall investigate the electronic properties of selected 
configurations. In particular, we first focus on the interplay of the magnetic behavior of 
the nh zig zag edges and of the spCC. We will then move on to describe the DFT-LSDA 
band structure, the vibrational properties, and the high-temperature stability of nh-C n 
configurations. 



4. Magnetism 



Zig-zag edges are generally known to be ferrimagnetic (Ml EE] due to non-bonding 
localized edge states. A detailed investigation of the magnetic properties of graphene 
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(a) nh-C 7 



(b) nh-C 8 



Figure 9. (Color online) The 0.01 hb/o,q magnetic isosurface of the ferromagnetic 
state of nh-C7 and of nh-Cg. 



edge in the context of a nanohole was carried out by Yu et al. [67]. That work focused 
on zig-zag edges (armchair ones are known to be nonmagnetic [69]), which made it 
convenient to study diamond- or hexagon-shaped holes with zig-zag edges only. The 
main conclusion of Ref. [67J regarding consecutive zig-zag edges is that the relative 
alignment of magnetic moments tends to be ferromagnetic when the edge atoms belong 
to the same graphene sublattice. This conclusion can be rephrased in terms of the angle 
between the two consecutive edges, which is defined as the angle between the in-plane 
outward vectors normal to the edges, as illustrated schematically Fig. [8] Ferromagnetic 
correlations occur when subsequent zig-zag edges are unrotated (0°) or rotated by 
120°, as would happen in a triangular hole. In the opposite case, the magnetization 
is antiferromagnetic, as occurs for zig-zag edges rotated by 60° or 180° (relevant, e.g. 
for a diamond or and hexagonal hole). 

Our rectangular nh involves two armchair edges, which are long enough to isolate 



rather effectively the magnetic moments localized at the two zig-zag edges. In Appendix 
|A| we study the magnetic properties of the edge of this nh. Following spCC insertion, 
all structures described in Sect. [3] preserve a nonzero absolute magnetization, associated 
to unpaired-spin electrons localized at the zig-zag nh edges. A significant magnetization 
is shared by the C n spCCs with odd n, while the even-n spCCs are non-magnetic, as 
illustrated by Fig. [9] for the nh-Cy and nh-Cs structures. 

The ferromagnetic structures of Fig. [9] are induced by the choice of a uniform 
starting magnetization used to initialize the electronic self-consistent calculation. To 
investigate other possible magnetic arrangements, we need to start off the self-consistent 
calculation with different magnetic arrangements of the individual atoms. To do this, 
we define two fictitiously different atomic species, both with the same chemical nature 
of C, but with initial magnetizations of opposite sign (±1 Bohr magneton). We place 
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Figure 10. (Color online) To investigate nonferromagnetic structures we consider 
two fictitious atomic carbon species fixing the spin polarization at the initial stage 
of the self-consistent calculation. C-j- (dark/red) and (clear/gray) carry positive 
and negative initial magnetization, respectively. Arcs mark all symmetry-independent 
nearest-neighbor Ising-type magnetic couplings Jy, see Eq. (Tip . 



these initially magnetically polarized atoms along the zig-zag edges in order to trigger 
the desired magnetic structure. Figure [10] illustrates one of many possible arrangements 
of the C-f (dark/red) and C± (clear/gray) atoms to initiate the self-consistent electronic- 
structure calculation. 

4-1- C 6 -nh 

We perform several self-consistent calculations for the nh-C6 structure, considering 



different starting magnetizations, as shown in Fig. 11, and determine the ground 



magnetic configuration. In agreement with Ref. [57], the ground-state configuration, 



Fig. |ll(a) has atoms of the the same magnetization in the same graphene sublattice 



(e.g. atoms labeled B, C, H, E in Fig. 10), and magnetization changes sign in passing 
from one sublattice to the other. The edge atoms bonded to even-n C n spCCs show 
little magnetism, mainly induced by the ferromagnetic interaction with neighboring 
atoms along the same zig-zag edge. 

This ground-state magnetic configuration is relaxed completely, and the resulting 
total energy E gs is taken as reference. Keeping fixed this fully relaxed ground atomic 
configuration, we repeat single self-consistent DFT-LSDA evaluations of the total energy, 
Stat, integrated magnetization M tot = f M z {r) d 3 f, and integrated absolute value of 
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magnetization M a b s = / |M 2 (r)| d 3 r, which we report next to each structure in Fig. 
The resulting individual magnetic configurations, violating the opposite-sublattice rule, 
are low-lying excitations, which we obtain in the DFT-LSDA simulations by changing 
appropriately the initial magnetizations of selected atoms. Relaxation of one of these 
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Figure 11. (Color online) Magnetization-density isosurfaces at +0.01/iB/ao 
(dark/red) and — O-OI/ib/oq (clear/gray), for several different nh-Cg magnetic 
structures. Frame (a): the DFT-LSDA ground state. Other frames: magnetically 
excited states. For each frame we report the relevant excitation energy, total 
magnetization and integrated absolute magnetization in Bohr magnetons /j,b- 
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Ising Parameter 


Value 


meV] 


Standard deviation 


meV] 


Eq — E gs 


2< 


SI 


5 


Jl 


-259 


5 


■h 




-4 


10 


Js 


-17 


10 



Table 2. (Color online) The interaction parameters of the Ising Model, Eq. ([lj, fitted 
on the DFT-LSDA values of the total energy of the magnetic configurations of nh-Cg 
zig listed in Fig. [XT] 



configurations shows very small displacements, not larger than 6 pm. 

The excitation energies of such states can be described approximately within a 
Ising-model scheme. The z component Si of the spin degree of freedom accounting for 
the magnetization at site i interacts with neighboring spins Sj, with an energy 

H S pi n = ~ ^ JijSiSj . (1) 



<i,j> 



According to the values of the absolute magnetization reported in Fig. 11 (a)[ it is 



appropriate to assume that each edge atom carries one Bohr magneton, i.e. one unpaired 
spin 1/2, thus Si = ±1/2. Accordingly, it makes sense to fit Ising-model energies only 
to configurations with an absolute magnetization significantly close to 8/xb- To avoid 
parameter proliferation, we neglect interactions between non-nearest-neighbor magnetic 
atoms. 



Figure 10 identifies the 3 independent Ising interaction parameters allowed by 
symmetry: J\ for the interactions between unpaired spins in different sublattices on 
edges rotated by 60° {J\ = Jab = Jcd = Jef = Jgh)] -h for the interactions within 
the same zig-zag edge, but "isolated" by the spCC {J% = Jbc = Jfg)\ an d J3 for the 
interactions across the armchair edge, representing opposite sublattices, or edges rotated 
by 180° (J3 = J ah = Jde)- We write the energy of a configuration as the sum of the 
magnetic energy -E sp in, approximated by the Ising expression ([!]), plus Eq, including 
covalency and all other interactions establishing the mean value of the total energy, 
averaged over all possible spin orientations. Specifically, for the nh-C6 zig structure, we 
have: 

-^spin — — J\ {SaSb + ScSn + SeSf + SqSh) (2) 
— J2 (SbSq + SpSc) — J3 (SaSh + SdSe) , 



so that, given the ground configuration of Fig. 11(a), we have E gs = E + E spin = 
E + J! + (J 3 - J 2 )/2. 

We estimate the Ising-model parameters by means of a linear fit of the DFT 
energies with expression ([2]). Table [2] reports the best-fit values of Eq and J^: all 
the exchange energies Jk turn out negative, reflecting antiferromagnetic interactions. 
The most significant value is J\, reflecting the strong antiferromagnetic coupling of 
adjacent unpaired spins on edge atoms belonging to opposite sublattices. J\ is over 
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Figure 12. (Color online) Comparison of the DFT-LSDA energy levels of the 



magnetic structures of Fig. 11 (left) with the spectrum (right) obtained using the 
Ising model, based on parameteres fitting the DFT values. 



one order of magnitude greater than the weakly antiferromagnetic coupling J3 across 
an armchair edge section. Given the fit standard deviation, the small value of J2 is 
compatible with null coupling. The obtained small negative value is the result of a 



strong cancellation between the energy-order reversed configurations of panels ll(c 



11(d), and the regularly ordered states of panels 11(a), 11(b), the latter matching the 
ordering Ref. [67] as expected. Indeed the Ising-model ground state and first-excited 
level turn out almost degenerate and actually in reversed order due to the small positive 



value of J2, as illustrated in Fig. 12 



In Fig. 12 , we compare the energies of the different configurations of Fig. 11 obtained 
by DFT calculation with those obtained using the fitted Ising model. A remarkable 
feature of the DFT excitation spectrum is the tiny splitting of the levels of panels 
which is hardly compatible with a simple nearest-neighbor Ising model. 



11 te- 



ll g 



Eventually Fig. [12] shows that the simple Ising model fails to describe the fine structure 
of the magnetic excitation of the edge atoms in the considered geometry. Only the 
significant J\ energy, fixing the rough structure of the spectrum, is determined with fair 
accuracy. Of course one could easily modify the model to include e.g. second-neighbor 
interactions, to fit the detailed level structure, but that would take any predictive power 
out of the model. 



4.2. C 7 -nh 

One may attempt a similar analysis for the odd-n spCCs, e.g. nh-Cy. As Fig. |9 (a) | shows, 
odd-n spCCs are magnetic, thus quite different from the even-n ones. This leads to two 
consequences for odd spCCs attached to nh: first the spin values at different sites are 
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different, and second the number of spin interactions to be considered is greater. This 
would leave little significance to a Ising model description. 

It is possible to at least identify the ground magnetic configuration, like we did 
for even-n chains. Figure [13] shows two different magnetic configuration of the nh- 



C7 straight structure. The ground-state configuration is the one of Fig. 13(a) , which 
follows the edge rules of Ref. [67]. The coupling between the spCC and nh edges is 
antiferromagnetic: this can be seen as a special case of the edge rules if the spCC atoms 
are seen as graphene atoms belonging to the edge but in the other sublattice relative 
to the outer zig-zag edge magnetic atoms. This coupling is so strong that it prevails 
over the weak antiferro long-range J3-type coupling. One can estimate this magnetic 
coupling energy between the end spCC atom and one of the nearest zig-zag edge atoms to 
approximately w 100 meV. The intra-spCC interaction is distinctly antiferromagnetic. 

5. Electronic Properties 

The magnetic properties described in the previous section are determined by the 
electronic structure. Before analyzing the nanohole-spCC system, it is useful to examine 
the simpler bands of a empty nh. We will track the the bands along the Brillouin- 



zone path shown in Fig. 14(a) . We sample the k-space path with points separated by 
1.75 x 10~ 13 nr 1 . 



Figure 14(b) displays the band structure near the Fermi energy for the empty 
nh superlattice of Fig. [4] One can identify two kinds of bands, with different spatial 
localization properties of their wave functions: (i) States localized at the hole edge 



(HE), such as the one depicted in Figs. 15(a). We use magenta crosses to track these 
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Figure 14. (Color online) (a): The Brillouin-zone (dashed) with the r — X — M — 
r — Y — M k-point path adopted for all band-structure calculations of the present 
paper, (b) Spin-majority Kohn-Sham band structure of the relaxed nh superlattice of 
Fig. |4j in the ferromagnetic configuration of Fig. 2(f) magenta crosses stand for HE 
bands localized on the hole-edge; green squares represent delocalized BU states. The 
plot focuses a 4 eV-wide energy region around the Fermi level (red dashed) for better 
readability. 





(a) A HE state 



(b) A BU state 



Figure 15. (Color online) Examples of k = electronic states of the nh. (a): a 
localized HE state; (b): a BU state. In Fig. 14(b) these states are located at the 
T-point crossing at —0.81 eV and 1.36 eV, respectively. 
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HE states in the bands-structure plots, such as Fig. 14(b) (ii) States like the one in 



Fig. 15(b) localized primarily on the bulk graphene atoms, with a weak component on 
the edge atoms. We label these states as BU, and identify them with green squares in 
band-structure plots. 

HE bands are generally flat, with little dispersion. The small but nonzero 
bandwidth of HE states is due to the residual interaction between the nh and its periodic 
images. A HE band touches the Fermi energy near X, and is therefore only partly filled, 
thus becoming the responsible of the edge magnetism discussed in Sect. |4j 

Liu et al. in their investigation of the band structures of a different graphene 
nanohole [70] , discovered the opening of band gaps for nanoholes with either armchair 
or zig-zag edges. In contrast, our graphene with a isolated nh exhibits no band gap at 
the Fermi energy, and retains the (semi) metallic character of graphene. Specifically, a 
BU metallic band crosses €f and shows a modest but distinct dispersion, with sizeable 
empty hole pockets near the X and Y points. 

Coming to the nh-C n systems, Figs. 16, 17, 18j and[l9]report details of the computed 
DFT-LSDA bands for the relaxed structures of Figs. [6] and [7j Also in these band 
structures we identify HE (magenta crosses) and BU (green squares) bands. In addition, 
band states significantly localized on the spCC atoms (CB) are identified by blue circles. 
One such state is depicted in Fig. 20(a) Occasional resonances of localized spCC and 



nh-edge states lead to hybrid localized states involving both, e.g. the one of Fig. |20(b) 



We identify such "CHE" states only in nh-Cg arm, and label them by black triangles 
in Fig. 18(a) Graphene bulk states often hybridize with the spCC molecular orbitals, 



thus acquiring a significant spCC components, as illustrated for example in Fig. 20(c) 



Our sample is too small to distinguish clearly between entirely localized states at the 
spCC/nh edge (whose bands would be perfectly flat in a realistically wide sample) and 
only partly localized hybrid states. 

Essentially all structures display a metallic behavior, due to one or several bulk 
bands crossing the Fermi energy. The different bonded spCCs affect the graphene nh 
bands quite considerably, by both shifting them and deforming them especially near 
the Fermi energy. In particular, the positions of several localized states at the nh edge 
change depending on the spCC state, and moreover spCC-specific localized states occur. 
For even-n spCCs, the CB states are energetically quite distant from the Fermi level, 
while odd-n spCCs exhibit a more metallic behavior, with spCC states quite close to the 
Fermi energy, and significant hybridization with the extended bulk states, consistently 
with results of Ref. [25l [7T] . 

The comparison of the spin majority and minority bands in Fig. 17 shows that 
magnetism affects the bulk bands only weakly. Magnetism appears to be associated to 
an energy shift of a few localized HE and (for odd spCCs) CB states near the Fermi 
level. The resulting effective exchange energy is ~ 0.3 eV. 

We perform several calculations of the nh-Cy structure for each of the considered 
spCC shapes: curved - Fig. 6(g) , s-curved - Fig. 6(h) , and straight - Fig. 6(i) All these 



geometries show basically identical band structures, e.g. the one reported in Fig. |18(b) 
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Figure 16. (Color online) Spin-majority band structures of (a) nh-Cs, represented 
in Fig. [6(a)] and of (b) nh-C 5 16, Fig. [6(b) 



Likewise, no special effect of the spCC curvature is apparent in the bands of nh-C 8 , 

is 



Fig. 19(a) Finally, the congestion of the bands near the Fermi level in Fig. 19(b) 
a consequence of the larger cell, and greater number of atoms and of electrons of this 
specific configuration. The general considerations (even-n spCC bands away from the 
Fermi energy, magnetism related to HE bands near the Fermi energy) apply also in this 
more intricate configuration. 
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(b) nh-Cg zig minority spin 



Figure 17. (Color online) Band structure for the majority (a) and minority (b) spin 



components of the ferromagnetic state of nh-C6 zig, depicted in Fig. 6(c 



6. Vibrational spectra 

We perform phonon calculation for a few stable structures of Sect. [3j We evaluate the 
phonon frequencies and eigenvectors of using standard density-functional perturbation 
theory, as implemented in the Quantum Espresso code p3l E2j- For comparison, the 
theoretical C-C stretching modes of polyynes C n H 2 (n = 8 — 12) evaluated with the 
same method match the experimental frequencies [73] to within 40 cm -1 . The size of 
the system is too large to evaluate the full dynamical matrix: although in principle 
possible, it would require a huge investment of computer time. 

We focus specifically on spCC "optical" stretching modes, which are prominent 
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Figure 18. (Color online) Spin-majority band structures of (a) the nh-Cg arm 
nanostructure, depicted in Fig. 6(d)| and (b) nh-Cy, Fig. 6(c) 



and characteristic in the experimental spectra of sp — sp 2 carbon in the spectral region 
near 2000 cm -1 , while all other vibrations (the "acoustic" spCC stretching modes, all 
bending modes, all graphene vibrations) overlap and lump together in a continuum 
extending from to 1600 cm -1 [231 [22] ■ We verified that the spCC stretching modes 
are influenced very little by faraway ligand atoms [73]. Accordingly, we only compute 
and diagonalize the part of the dynamical matrix at T relative to displacements of the 
atoms of the spCC, plus its first and second neighbors in the graphene sheet. The error 
in the vibrational frequency induced by this approximation can be estimated ~ 1 cm -1 . 

By analyzing the displacement pattern of the normal modes of the even-n spCCs, 
it is straightforward to identify the "a" modes characterized by the strongest Raman 
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(b) nh-2C 6 



Figure 19. (Color online) Spin-majority band structures of nh-Cs, see Fig. 6(f) and 
of nh-2C 6 , see Fig. [7] 



Structure name 


Raman frequencies 


cm 1 




IR frequencies [cm 1 




nh-C 5 
nh-C 6 zig 


1323, 1368 
1777, 1878 


1332 
1939 



Table 3. Wavenumber of Raman and IR spCC stretching frequencies calculated for 
the nh-C n structures. The most intense Raman and IR frequencies are highlighted in 
bold. 
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A CB state 



(b) A CHE state 



(c) A resonant bulk state 



Figure 20. (Color online) Like in Fig. 15 but for nh-Cg zig and nh-Cg arm. (a): 
a typical CB state localized mainly on the polyyne - the T state at energy —1.16 eV 
in Fig. 17(a)| (b): a CHE state, localized jointly on the polyyne and the nh edge - 
0.50 eV in Fig. [18(a) 



the T state at energy 
extension on the polyyne 



(c): a bulk state showing a significant 



the T state at energy —0.46 eV in Fig. 17(a)| 



and IR absorption [251 EH EEl [75] . We take advantage of this pattern recognition to 
avoid a computationally expensive explicit evaluation of the Raman and IR intensities. 
Table [3] reports the computed frequencies of the spCC stretching modes, with the 
most intense Raman and IR mode highlighted. Note that the wavenumbers of the 
frequencies are significantly lower than the characteristic stretching-frequencies of free 
spCCs (1950 — 2300 cm -1 ). The reason is the tensile strain to which the spCCs are 
subjected by binding to the nh edges. In nh-Cs, the length of the chain, including the 
bonds between the spCC and the nanohole, is 15% longer than the isolated spCC; this 
elongation leads to frequencies much softer than typical polyyne ones. The chain length 
in nh-C6 zig is only 5% longer than isolated length, and the frequencies come much 
closer to the typical frequencies of free spCCs. 

The results of the present section do not imply that spCCs in a context of 
nanostructured sp-sp 2 carbon material should vibrate at much different frequencies 
from their molecular counterparts [54J . Quite on the contrary, previous calculations and 
experiments confirm that fully relaxed spCCs terminated by sp 2 -type material exhibit 
very similar frequencies to those of molecular spCCs [251 EH E5]- The results of the 
present calculations suggest instead that unrelaxed tensile strain in nanostructured sp 
- sp 2 carbon material is likely to induce significant frequency shifts of the spCC modes. 
Depending on the method of production of spCC-containing material (e.g. cluster beam 
formation/deposition [23] vs. atomic wires stretched out from pulled graphene sheets 
[27J[2H]), whenever a sizeable tensile strains remain frozen in the sample, one is to observe 
a corresponding distribution of the observed vibrational frequencies, quite independent 
of the frequency shifts associated to different lengths and terminations of the spCCs 
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7. High-temperature stability 

All studied configurations are local minima of the adiabatic potential energy, i.e. 
metastable allotropes of carbon. Given sufficiently long time, the spCCs are expected to 
degrade, for example by recombining with the nh edge and extend energetically favored 
sp 2 graphene. This possibility is however very remote at low temperature, because in 
this interconversion process the energy barriers to be crossed are quite substantial. 

Since several detailed types of sp — > sp 2 processes are possible, a full study of the 
spCC degradation is beyond the scope of the present paper. We content ourselves with 
a semi-quantitative estimate of the thermal stability and the degradation mechanisms 
of spCCs bonded to carbon sp 2 nanostructures by running comparably long high- 
temperature molecular-dynamics (MD) simulations, and monitor the eventuality of 
spCC decomposition as a function of the simulated temperature. 

As, due to the size dependency of statistical fluctuations, the longer the spCCs the 
higher is the chance of chain breaking. We therefore prefer to simulate a larger version 
of the model of Sect. |3j namely a Cio chain bound to a nh large enough for it to fit 
loosely. Also the graphene plane is represented by four, rather than three hexagonal rings 
separating the nh periodic replicas. The resulting nh-Cio structure involves 184 atoms 
and 736 electrons. Even using Car-Parrinello dynamics, it would be a formidable task 
to simulate several sufficiently long runs to expect a significant chance to observe spCC 
dissociation at a temperature comparable to experiment. For the present task therefore 
we abandon the ab-initio DFT-LSDA method for the treatment of the electronic degrees 
of freedom, and replace it with a tight-binding (TB) model [76]. We adopt the TB 
scheme of Xu et al. [77], which has been applied successfully to investigate several low- 
dimensional carbon systems [78], [79] EH El] , as it reproduces well the experimental bulk 
equilibrium distance G? grap h = 1.4224 A of graphite, and the bulk structure and elastic 
properties of sp 3 diamond. Since all interatomic interactions vanish at a cutoff distance 
r c = 2.6 A, which is shorter than the interlayer spacing of graphite, 3.35 A [82], we focus 
on a single-layer model like we did for the DFT-LDA model above. 

To validate the TB force field for spCCs, we compute the total energy of a Cu 
linear chain as a function of the length <i 7 _ 8 of its central bond (between atoms 7 and 
8), which we keep fixed while all other bonds are allowed to relax. Figure 21 compares 
this total energy, referred to its value at full relaxation, as obtained with DFT-LSDA 
and with TB. The two models exhibit significant differences, in particular the TB model 
has a level crossing to a dissociative regime above 1.9 A, while nothing of the sort occurs 
in the DFT-LDA band-structure calculation, which takes care of level degeneracies by 
selecting a spin-1 magnetic state. This problem with the TB model is characteristic 
of unsaturated conditions, while, in close-shell electronic configurations, dissociation is 
more regular. In general, the TB model is likely to be more "fragile" than a more 
realistic DFT-LDA. In turn the latter is also expected to be less strongly bounded than 
real spCCs, due to missing long-range attractive Van-der-Waals polarization correlation 
effects. We must then conclude that all quantitative stability evaluations based on the 
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Figure 21. (Color online) Comparison of the total energy variation for a C14 linear 
chain as a function of the (constrained) length dj-% of its central bond (between atoms 
7 and 8), as obtained with DFT-LSDA and with TB. 



TB model are underestimates of the actual stability in experiment. Periodic boundary 
conditions (matching the zero-temperature lattice parameter of graphene) suppress long 
wavelength fluctuations, resulting in a stabilization of the thermal fluctuations, which 
in the thermodynamical limit would make the ID - 2D structure unstable. 

We run microcanonical (constant-energy) TB molecular dynamics (TBMD) 
simulations, since the temperature fluctuations are small enough (~ 10%) in a sample of 
this size for temperature to be considered a fairly well defined quantity. The advantage 
of the microcanonical ensemble is that no thermostat artifacts, and in particular no 
dissipative term as in Langevin or Nose-Hoover thermostats, can affect the atomic 
dynamics, allowing for a full account of local fluctuations to produce whatever bond 
breaking they may lead to. A disadvantage of the constant-energy MD is that, if a 
significant bonding breakdown of a part of the nanostructure occurs, the corresponding 
potential-energy increase occurs at the expense of the kinetic energy, thus the system 
may artificially cool down, thus hindering further decomposition. In practice, this 
problem has little importance for the system size considered. We use a time step of 
0.5 fs, small enough to guarantee a rigorous global energy conservation within 0.02 eV, 
or 0.001%. 

We run extended simulations of the model nh-Cio structure at different 
temperatures, starting with a sampling of initial conditions. We generate starting 
configurations by beginning with the fully relaxed configuration and running three 
successive equilibration runs (0.1 ps, 0.1 ps and 0.5 ps), with randomized initial 
velocities, taken from a Gaussian distribution matching the Boltzmann distribution 
at the target temperature. In the figures we indicate the resulting temperature, 
with an error bar appropriate to the determination of the average along the whole 
simulation, (thus not estimating the instantaneous temperature fluctuations, which are 
much larger). In the determination of this average temperature, we drop the first 100 fs, 
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(a) t = ps 



(b) t = 1.5 ps 



(c) t = 2.5 ps 



(d) t = 2.7 ps 



Figure 22. (Color online) Successive snapshots of a sample TBMD simulation at a 
temperature T = (3991 ± 10) K of nh-Ci . 



---- 3921 KB 
- 4006 K 




Figure 23. (Color online) The time dependence of the longest bond length d n 
Eq. ([3j), of the Cio chain in four simulations at different temperature. 



to minimize the systematic oscillations induced by starting with initial random velocities 
with little correlation to the forces. 



Figure [22] displays successive frames of an example simulation illustrating that edge 

often occur before the earliest spCC breakdown 



reconstruction processes 
event 



Fig. 22(b) 



Fig. 22(d) To monitor these breakdown events, a clear indicator is the longest 
C-C bond length relative to the 11 bonds of the Cio chain and of the chain ends to the 
attached nh edge atoms: 



(L 



max dn—c-A 
=0,1, ..io 1 H 



(3) 



Figure [23] reports the time dependency of d max following 4 independent simulations 
carried out at different temperature, with different initial states. Chain breakdown, as 
happens between frames 22(c)| and 22(d)| is signaled by the rapid increase of one of the 
bond lengths beyond 5 A. The spCC breakdown may be followed by recombination of 
the chain into the nh edge, or even expulsion of a section of the spCC into vacuum. 
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Figure 24. (Color online) Circles: the average inverse time before the breakdown 
of the Cio spCC, as a function of inverse temperature. Dashed line: a Arrhenius fit, 
Eq. @. 



As suggested by Fig. 23 spCC breakdown occurs, on average, earlier and earlier 
for increasing temperature. By repeating the numerical simulations for different initial 
conditions but similar temperature we estimate an average decay rate by averaging the 
inverse times before decay. The average decomposition rate time k is an increasing 
function of temperature. In simulations done at substantially lower temperature than 
the ones considered in Fig. 23 one would need to wait too long to observe decomposition, 
while at much higher temperature decomposition occurs immediately after start. 

If one can assume that one type of process (bond breaking) dominates over all 
decomposition channels, the decomposition rate k is expected to be in the Arrhenius 
form 

* = (4) 

We estimate the attempt frequency A and the effective energy barrier E a of the TBMD 
model by fitting the Arrhenius plot in Fig. 24 The value E a = 4.2 eV matches the 
TB breakup of Fig. 21 this effective activation barrier is surely an underestimation of 
the actual barrier against breakup. The estimated attempt rate A = 1.5 x 10 17 s -1 , 
although probably slightly overestimated, reflects the high number of breakup channels 
available for decomposition of the spCC [83J. If we assume the computed values for 
E a and A, we extrapolate thermal decay rates of spCCs in sp 2 carbon of the order of 
k = 7.5 x 1CT 5 s" 1 at 1000 K, and k = 1.5 x 10" 54 s" 1 at room-temperature (300 K). 
Overall, the calculations of the present section confirm a substantial stability of spCCs 
in the solid state and in vacuum. In the lab, whenever sp — sp 2 carbon is not kept in 
vacuum, chemical decay mechanisms are therefore likely to overcome the thermal ones. 
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8. Discussion and conclusion 



The present work collects extensive investigation of the geometry, electronic structure, 
magnetic properties, and dynamical stability of spCCs attached to sp 2 graphitic 
fragments. When a spCC binds to zig-zag graphene edges, its polyynic character is 
attenuated to a value intermediate between those typical of cumulenes and polyynes. 
The attachment of a spCC to the graphene edge is very stable: we predict stabilization 
energies near 6 eV per bond between each chain end and the sp 2 regions. Thermal 
excitations typically break bonds along the spCC with similar probability to those 
formed with the graphene edge, which indicates a very solid attachment. 

Odd spCCs in the nh display a metallic behavior, with at least one band pinned 
to the Fermi energy, while even spCCs have little overlap with the states at the Fermi 
level [501 ED EH EH]- The partly filled states of odd-n spCCs are associated to nonzero 
magnetization related to a spin triplet state of the n bonds. Even-n spCCs are instead 
insulating and non-magnetic, and in this context only the magnetic moments of the 
graphene edge contribute to the magnetism of the nh-C2 m structures, and only the bulk 
graphene states provide conducting bands. 

We compute also the vibrational modes, and specifically the optical C=C stretching 
modes which emerge as a characteristic signature of spCCs in Raman and IR 
spectroscopies. Our calculations show that the vibrational frequencies can be quite 
substantially red-shifted when spCCs are kept under tensile stress. Indeed, the weaker 
stability and correspondingly faster decomposition rate of strained spCCs is likely to 
play a significant role in the overall blue shift of the sp-carbon peak in the Raman 
spectrum of the decaying spCCs in cluster-assembled sp-sp 2 film [25J. 
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Appendix A. Magnetism at the edge of the empty nh 

As magnetism is intrinsic of the zig-zag graphene edge, thus even a empty nanohole 
exhibits a range of magnetic states similar to those arising in the presence of even-n 



spCCs. The zig-zag edge atoms involved in magnetism are clearly identified in Fig. Al 



Also here, each atomic site A-L carries a magnetization close to 1 ub- As illustrated in 



Fig. |A1[ even though the number of spin-carrying atoms is larger than in the case of 



Sect. 4.1, only three independent nearest-neighbor interactions need to be considered. 



Two of them, J\ and J3, represent interaction between two edges belonging to different 
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Figure Al. (Color online) Magnetization-density isosurfaces at +0.01 \ib /«o 
(dark/red) and —0.01 ^b/cIq (clear gray) for the magnetic ground state of the unpaired- 
spin electrons localized at the zig-zag edges of the empty nh. Letters label the spin- 
carrying atomic sites. Arcs mark all symmetry-independent nearest-neighbor Ising- 
type magnetic couplings Ji, see Eq. ([!]). 



Ising Parameter 


Value 


meV] 


Standard deviation 


meV] 


Eq — E gs 


300 


12 


Jl 


-248 


12 


■h 


10 


12 




-53 


19 



Table Al. (Color online) The individual Ising-model parameters computed for the 
nh magnetic configurations. 



sublattices (where we expect an anti-ferromagnetic character). In contrast, the coupling 
J2 accounts for the interaction of spins within the same edge, and is therefore expected 
to be ferromagnetic, according to Ref. [57] . 

As done in Sect. |4.1| for the nh-Cg structures, we make a linear fit of all considered 
magnetic configurations (shown in Fig. A2), to evaluate the values of the Ji and of Eq. 



In the Ising model, the total energy is written as: 

Etot = E + E spin = Eq — Ji {SaSb + SdSe + SpSc + SkSl) + 

— J 2 (SbSq + ScSd + SqSh + ShSk) + 
-J 3 (S A S L + S E S F ) . (A.l) 

The result of the linear fit of the energies of the magnetic configurations of Fig. |A2| is 



reported in Table Al Like in nh-C6, J\ is much larger than the other antiferromagnetic 
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(a) Ground State 

E tot = E gs = -13017.539 eV 

M tot = 0.00 

M abs = 11.76 




(d) E to t — E gs 
Mtot = 2.00 
M abs = 10.77 



86 meV 



(b) E^t 
M^t = 0.00 
M abs = 10.05 



E gs + 52 meV 



(c) E toi — Egs 
M^t = -2.00 
M abs = 9.71 



72 meV 





(e) E tot = E gs + 513 meV 
Mtot = 0.00 
M abs = 11.52 



(f ) -^tot — J^gs 



= E gs + 565 meV 
M^t = 10.00 
M abs = 10.77 



Figure A2. 

(clear/gray) 



(Color online) Magnetic positive (dark/red) +0.01/ib/oo an£ l negative 
-0.01/ib/oo isosurfaces for the empty nanohole structures. 



coupling J 3 . Indeed, the J\ interaction is very similar to the one obtained in the 
calculations with the nh-C6 structure, see Table |2j Moreover, like in the nh-C6 case, 
given the relevant error bar the weakly ferromagnetic J<i is in fact compatible with a 
null value, which is somewhat surprising for neighboring atoms along the same zig-zag 
edge. 



Figure A3 compares the DFT-LSDA energy levels and those obtained using the 
fitted Ising model. Here the fit agrees better than in the nh-C6 case, but clearly the role 
of the J 2 coupling is contradictory, which explains its small value. 
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